
library(survival)
clin<-read.delim("gbm_survival_6_21_2011_parsed" ,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="",
header=TRUE,
col.names=c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE"),as.is=c(TRUE,TRUE,FALSE,TRUE,TRUE,TRUE,TRUE,TRUE,TRUE,FALSE),colClasses=c("character","numeric","factor","numeric","numeric","numeric","integer","integer","integer","factor"))
#!!! Note Edited data from BORA
#e <- read.delim("merge_expr",header=TRUE,sep="\t",quote="\"",dec=".",fill=TRUE,comment.char="",na.strings="NA")
e <- read.delim("combined_GBM_normal_final_maf05_manhattan.tped",header=TRUE,sep=" ",quote="\"",dec=".",fill=TRUE,comment.char="",na.strings="NA")
esamples <- colnames(e)
geneids <- rownames(e)
# rownames do not tolerate "-", so the "-" were auto-converted to "."
clin$SUBJECT<-gsub("-",".",clin$SUBJECT)


res<-data.frame(gene=c(rep("",length(geneids))),pgene=c(rep(0.0,length(geneids))),pmodel=c(rep(0.0,length(geneids))),effect=c(rep(0.0,length(geneids))),stringsAsFactors=FALSE)
rownames(res)<-geneids
# Note: First, fit 1 clinical term at a time to decide wether to include in the model.
#	fit1 <- coxph(Surv(OS,ALIVEOS) ~ EXPR+AGE+GENDER+SITE,data=dt)


do.test <- function(i) {
   g<-geneids[i]
   et<-t(e[g,])
   dt <- merge(clin,et,by.y="row.names",by.x="SUBJECT")
   colnames(dt) <- c("SUBJECT","AGE","GENDER","OS","PFS","DFS","ALIVEOS","ALIVEPFS","ALIVEDFS","SITE","EXPR")
   fit1 <- coxph(Surv(OS,ALIVEOS) ~ EXPR+AGE,data=dt)
   s<-summary(fit1)$coefficients
   e_signif<-as.numeric(s["EXPR","Pr(>|z|)"])
   effect<-as.numeric(s["EXPR","coef"])
   model_signif<-as.numeric(summary(fit1)$logtest["pvalue"])
   return(c(g,e_signif,model_signif,effect))
}
res<- do.call(rbind, lapply(seq_along(geneids), FUN=do.test))
colnames(res) <- c("gene","pgene","pmodel","effect")
write.table(res,file="genotyping_survival.txt",sep="\t",quote=FALSE,row.names=FALSE)

